#!/usr/bin/env julia
"""
Manufactured Solution Validation Example

This example validates the NSEMSolver.jl implementation against manufactured
solutions to verify code correctness and convergence behavior.
"""

using NSEMSolver

function validate_manufactured_solution_2d()
    println("🧪 Manufactured Solution Validation (2D)")
    println("=" * 50)
    
    # Create manufactured solution for testing
    manufactured_sol = create_manufactured_solution(2)
    
    # Validate the manufactured solution properties
    println("📋 Validating manufactured solution properties...")
    is_valid = validate_manufactured_solution(manufactured_sol)
    
    if is_valid
        println("✅ Manufactured solution validation passed!")
    else
        println("❌ Manufactured solution validation failed!")
        return false
    end
    
    # Test convergence for multiple grid refinements
    println("\n🔍 Testing grid convergence...")
    
    N_values = [3, 4, 5]  # Polynomial orders to test
    errors = Float64[]
    
    for N in N_values
        println("\n  Testing N = $N...")
        
        options = NSOptions(
            N = N,
            n_block = 2,             # Small domain for quick testing
            nu = 0.01,
            tfinal = 0.1,            # Short time for validation
            cfl = 0.3,
            tol = 1e-8,
            solver = :julia,
            method = :sem,
            adaptive_refinement = false,  # Disable for convergence test
            verbose = false           # Quiet for comparison
        )
        
        # Solve the system
        result = solve_navier_stokes_2d(options)
        
        if result.converged
            # Calculate error against manufactured solution
            error = compute_manufactured_solution_error(
                manufactured_sol,
                result.u, result.v, nothing, result.p,
                result.multidomain.xx, result.multidomain.yy, nothing,
                options.tfinal
            )
            
            total_error = error.total_l2_error
            push!(errors, total_error)
            
            println("    L2 error: $(@sprintf("%.2e", total_error))")
            println("    u L2 error: $(@sprintf("%.2e", error.u_l2_error))")
            println("    v L2 error: $(@sprintf("%.2e", error.v_l2_error))")
            println("    p L2 error: $(@sprintf("%.2e", error.p_l2_error))")
        else
            println("    ❌ Solution did not converge for N = $N")
            push!(errors, Inf)
        end
    end
    
    # Analyze convergence rate
    if length(errors) >= 2 && all(errors .< Inf)
        println("\n📊 Convergence Analysis:")
        for i in 2:length(errors)
            if errors[i] > 0 && errors[i-1] > 0
                rate = log(errors[i-1] / errors[i]) / log((N_values[i] + 1) / (N_values[i-1] + 1))
                println("  N=$(N_values[i-1]) → N=$(N_values[i]): convergence rate ≈ $(@sprintf("%.2f", rate))")
            end
        end
        
        if errors[end] < 1e-4
            println("✅ Grid convergence test passed!")
            return true
        else
            println("⚠️ Final error may be higher than expected")
            return false
        end
    else
        println("❌ Grid convergence test failed!")
        return false
    end
end

function test_manufactured_solution_properties()
    println("\n🔬 Testing Manufactured Solution Properties")
    println("-" * 40)
    
    # Test both 2D and 3D manufactured solutions
    for dim in [2, 3]
        println("\n  Testing $(dim)D manufactured solution...")
        
        sol = create_manufactured_solution(dim)
        
        # Test evaluation at specific points
        test_points = dim == 2 ? [(0.5, 0.5, 0.0)] : [(0.3, 0.7, 0.2)]
        test_time = 0.5
        test_nu = 0.01
        
        for point in test_points
            x_val, y_val, z_val = point
            
            # Test velocity evaluation
            u_val = calculate_manufactured_u(sol, x_val, y_val, z_val, test_time)
            v_val = calculate_manufactured_v(sol, x_val, y_val, z_val, test_time)
            w_val = calculate_manufactured_w(sol, x_val, y_val, z_val, test_time)
            p_val = calculate_manufactured_p(sol, x_val, y_val, z_val, test_time)
            
            # Test forcing evaluation  
            fx, fy, fz = calculate_manufactured_forcing(sol, x_val, y_val, z_val, test_time, test_nu)
            
            println("    At ($x_val, $y_val, $z_val), t=$test_time:")
            println("      u = $(@sprintf("%.3f", u_val))")
            println("      v = $(@sprintf("%.3f", v_val))")
            if dim == 3
                println("      w = $(@sprintf("%.3f", w_val))")
            end
            println("      p = $(@sprintf("%.3f", p_val))")
            println("      fx = $(@sprintf("%.3f", fx))")
            println("      fy = $(@sprintf("%.3f", fy))")
            if dim == 3
                println("      fz = $(@sprintf("%.3f", fz))")
            end
        end
        
        # Test divergence-free condition
        is_div_free = validate_manufactured_solution(sol)
        println("    Divergence-free: $(is_div_free ? "✅" : "❌")")
    end
    
    return true
end

function benchmark_solvers()
    println("\n⚡ Solver Benchmark Comparison")
    println("-" * 40)
    
    # Create a medium-sized problem for benchmarking
    options_base = NSOptions(
        N = 4,
        n_block = 3,
        nu = 0.01,
        tfinal = 0.5,
        cfl = 0.4,
        tol = 1e-6,
        adaptive_refinement = false,
        verbose = false
    )
    
    solvers = [:julia]
    if HAS_PETSC
        push!(solvers, :petsc)
    end
    if HAS_GCR
        push!(solvers, :gcr)
    end
    
    results = Dict{Symbol, Any}()
    
    for solver in solvers
        println("\n  Benchmarking $(solver) solver...")
        
        options = NSOptions(options_base; solver = solver)
        
        start_time = time()
        result = solve_navier_stokes_2d(options)
        end_time = time()
        
        if result.converged
            results[solver] = (
                time = end_time - start_time,
                iterations = result.iterations,
                residual = result.residual_norm,
                success = true
            )
            
            println("    Time: $(@sprintf("%.3f", results[solver].time)) s")
            println("    Iterations: $(results[solver].iterations)")
            println("    Final residual: $(@sprintf("%.2e", results[solver].residual))")
        else
            results[solver] = (success = false)
            println("    ❌ Failed to converge")
        end
    end
    
    # Find fastest solver
    successful_solvers = [k for (k, v) in results if haskey(v, :success) && v.success]
    if length(successful_solvers) > 1
        fastest = argmin([results[k].time for k in successful_solvers])
        fastest_solver = successful_solvers[fastest]
        println("\n  🏆 Fastest solver: $(fastest_solver) ($(@sprintf("%.3f", results[fastest_solver].time)) s)")
    end
    
    return results
end

function main()
    println("🧪 NSEMSolver.jl - Manufactured Solution Validation")
    println("=" * 60)
    
    success = true
    
    # Test 1: Validate manufactured solution properties
    try
        test_manufactured_solution_properties()
        println("✅ Manufactured solution property tests passed")
    catch e
        println("❌ Manufactured solution property tests failed: $e")
        success = false
    end
    
    # Test 2: Validate grid convergence
    try
        convergence_passed = validate_manufactured_solution_2d()
        if !convergence_passed
            success = false
        end
    catch e
        println("❌ Grid convergence validation failed: $e")
        success = false
    end
    
    # Test 3: Benchmark solvers
    try
        benchmark_results = benchmark_solvers()
        println("✅ Solver benchmarking completed")
    catch e
        println("❌ Solver benchmarking failed: $e")
        success = false
    end
    
    println("\n" * "=" * 60)
    if success
        println("🎉 All validation tests passed!")
        println("✅ NSEMSolver.jl implementation is verified against manufactured solutions")
    else
        println("⚠️  Some validation tests failed!")
        println("❌ Please check the implementation for potential issues")
    end
    
    return success
end

# Run the validation
if abspath(PROGRAM_FILE) == @__FILE__
    main()
end